/*=+--+=#=+--     Unmanned Aerial System Management Software      --+=#=+--+=#*\
|          Copyright (C) 2011 Regents of the University of Colorado.           |
|                             All Rights Reserved.                             |

     This program is free software: you can redistribute it and/or modify
     it under the terms of the GNU General Public License version 2 as
     published by the Free Software Foundation.

     This program is distributed in the hope that it will be useful,
     but WITHOUT ANY WARRANTY; without even the implied warranty of
     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     GNU General Public License for more details.

     You should have received a copy of the GNU General Public License
     along with this program.  If not, see <http://www.gnu.org/licenses/>.

            Jack Elston                       Cory Dixon                        
|           elstonj@colorado.edu              dixonc@colorado.edu              |
|                                                                              |
\*=+--+=#=+--                 --+=#=+--+=#=+--                    --+=#=+--+=#*/

/***********************************************************************
 *
 * FILENAME:
 * gtk_gpsmap.c
 *
 * PURPOSE:
 *	Provide a GPS map interface to the user.
 *
 * LAST MODIFIED:
 *
 ***********************************************************************/

#include "geoImage.h"
#include "utilities.h"

#include <iomanip>
using namespace std;

GeoImage::GeoImage(void)
{
	tiff	= NULL;
	gtif	= NULL;
	pixbuf = NULL;

	width = height = 0;
	dx_dPp = dy_dPp = 0.0;
	img_tieX = img_tieY = m_tieX = m_tieY = m_resX = m_resY = 0.0;
}

GeoImage::~GeoImage(void)
{
	Close();
}

void GeoImage::Close(void)
{
	if (gtif) GTIFFree(gtif);
	gtif = NULL;
	if (tiff) XTIFFClose(tiff);
	tiff = NULL;

	if(pixbuf) gdk_pixbuf_unref(pixbuf); 
	pixbuf = NULL;

	//if(imageProj)
		//pj_free(imageProj);
}

bool GeoImage::Open(String file, String gfw )
{
	filename    = file;
	wf_filename = gfw;

	return Open();
}

bool GeoImage::Open()
{
	if( filename == "" )
		return false;

	if( gtif || tiff || pixbuf)
		Close();

	// test for valid image file and get info
	GdkPixbufFormat *fmt = gdk_pixbuf_get_file_info(filename.c_str(), &width, &height);
	if( fmt == NULL ) {
		cout << "GeoImage::Open() - unrecognized format" << endl;
		return false;
	}

	//cout << "Loading " << fmt->name << " image ...";
	cout.flush();

	// see if there is a world file. if so, just load the image
	// otherwise see if it is a geoTIFF
	if( readWorldFile() )
		return image_open();
	else if( filename.find(".tif") > 0  )
		return geoTIFF_open();

	cout << " did not find a GEO referenced image" << endl;
	return false;
}

bool GeoImage::image_open()
{
	GError *error = NULL;

	pixbuf = gdk_pixbuf_new_from_file (filename.c_str(), &error);
	if( error ) {
		printf( "GeoImage::image_open - ERROR: %s", error->message);
		return false;
	}

	// manualy calculate deg/pixel
	calcDegPerPixel();

	//printImageInfo();

	//if( fabs(dx_mPp) != fabs(dy_mPp) ) 
		//cout << " !!! This is not a square (rectified) image. It WILL NOT be displayed correctly" << endl;

	return true;
}

bool GeoImage::geoTIFF_open()
{
	// Open TIFF descriptor to read GeoTIFF tags 
	tiff = XTIFFOpen(filename.c_str(),"r");
	if (tiff)
	{
		// Open GTIF Key parser; keys will be read at this time
		gtif = GTIFNew(tiff);
		if (gtif)
		{
			//const int key_size = 2048;
			//char key_data[key_size];
			//memset(key_data,0,(size_t)key_size);

			double * d_list = 0;
			int d_list_count = 0;

			// make a  pixbuf copy, it might be dumb, but the easies way is to simply
			// reload the file again
			GError *error = NULL;
			pixbuf = gdk_pixbuf_new_from_file(filename.c_str(), &error);
			if( error ) {
				printf("GeoImage::geoTIFF_open - ERROR: %s\n", error->message);
				Close();
				return false;
			}

			// Get the GeoTIFF deffinition value
			if( !GTIFGetDefn( gtif, &defn ) ) {
				printf("GeoImage::geoTIFF_open - could not get geotiff definition structure");
				Close();
				return false;
			}

			if( defn.GCS != KvUserDefined )
			{
				char	*pszName = NULL;
				short pnDatum, pnPM, pnUOMAngle;
				GTIFGetGCSInfo( defn.GCS, &pszName, &pnDatum, &pnPM, &pnUOMAngle);
				if( pszName == NULL )
					pszName = CPLStrdup("(unknown)");

				gcsSTR = pszName;
				CPLFree( pszName );
			}

			switch( defn.Model ) {
				case ModelTypeProjected:   /* Projection Coordinate System */
					modelSTR = "PROJ";
					break;
				case ModelTypeGeographic:  /* Geographic latitude-longitude System */
					modelSTR = "GEO-LL";
					break;
				case ModelTypeGeocentric:  /* Geocentric (X,Y,Z) Coordinate System */
					modelSTR = "GEO-XYZ";
					break;
				default:
					modelSTR = "UNKNOWN";
			}

			if (TIFFGetField(tiff, GTIFF_TIEPOINTS, &d_list_count, &d_list))
			{
				img_tieX = d_list[0];
				img_tieY = d_list[1];

				m_tieX = d_list[3];
				m_tieY = d_list[4];
			}
			if (TIFFGetField(tiff, GTIFF_PIXELSCALE, &d_list_count, &d_list))
			{
				m_resX = d_list[0];
				m_resY = d_list[1];
			}

			// image width and height
			TIFFGetField( tiff, TIFFTAG_IMAGEWIDTH, &width );
			TIFFGetField( tiff, TIFFTAG_IMAGELENGTH, &height );

			// manualy calculate deg/pixel
			calcDegPerPixel();

			// dump the GeoTIFF metadata to std out 
			//GTIFPrint(gtif,0,0);
			//PrintCorners(true, false);
			//printImageInfo();

			//cout << "Model=" << defn.Model << endl
			     //<< "PCS=" << defn.PCS<< endl
			     //<< "GCS=" << defn.GCS  << endl
			     //<< "Datum=" << defn.Datum << endl
			     //<< "Ellipsoid=" << defn.Ellipsoid << endl
			     //<< "ProjCode=" << defn.ProjCode << endl
				  //<< "MapSys=" << defn.MapSys << endl
				  //<< "Zone=" << defn.Zone << endl;
	
			//if( fabs(dx_mPp) != fabs(dy_mPp) ) 
				//cout << " !!! This is not a square (rectified) image. It may not be displayed correctly" << endl;

			cout << endl;

		} else
			cout << "GeoImage::Open() - could not read as GEOTiff image: " << filename << endl;
	} else
		cout << "GeoImage::Open() - could not open as tiff image: " << filename << endl;

	return true;
}

bool GeoImage::calcDegPerPixel()
{
	if( gtif ) {

		if( defn.Model == ModelTypeGeographic ) {
			if( defn.UOMAngle == Angular_Degree ) {
				dx_dPp = m_resX;
				dy_dPp = m_resY;
			} else if ( defn.UOMAngle == Angular_Radian ) {
				dx_dPp = m_resX * 180.0/M_PI;
				dy_dPp = m_resY * 180.0/M_PI;
			} else {
				cout << "GeoImage::calcDegPerPixel() - can't convert all angular units yet" << endl;
				return false;
			}     

			// since coordinates are angular, calculate linear resultion
			dx_mPp = LON_TO_M( dx_dPp, m_tieY );
			dy_mPp = LAT_TO_M( dy_dPp );

		} else if( defn.Model == ModelTypeProjected ) {
			if( defn.UOMLength == Linear_Meter ) {
				dx_mPp = m_resX;
				dy_mPp = m_resY;
			} else {
				cout << "GeoImage::calcDegPerPixel() - can't convert all linear units yet" << endl;
				return false;
			}

			double x,y;	
			getTiePoints(x,y);

			// since coordinates are linear, calculate angular resolution
			dx_dPp = M_TO_LON( dx_mPp, y );
			dy_dPp = M_TO_LAT( dy_mPp );

		} else 
			return false;

	} else {
		dx_dPp = m_resX;
		dy_dPp = m_resY;

		dx_mPp = LON_TO_M( dx_dPp, m_tieY );
		dy_mPp = LAT_TO_M( dy_dPp );
	}

	//cout << " dx=" << dx_dPp << " [deg/pixel] " << " dy=" << dy_dPp << " [deg/pixel] " << endl
	     //<< " mx=" << dx_mPp << " [m/pixel] " << " my=" << dy_mPp << " [m/pixel] " << endl;

	return true;
}

/* 
 * Report the file(s) corner coordinates in projected coordinates, and 
 * if possible lat/long.
 */	
bool GeoImage::ReportACorner(const char * corner_name, double x, double y, int inv_flag, int dec_flag )
{									  
	double	 x_saved, y_saved;

	if( !gtif )
		return false;

	/* Try to transform the coordinate into PCS space */
	if( !GTIFImageToPCS( gtif, &x, &y ) )
		return false;							

	x_saved = x;
	y_saved = y;								

	printf( "%-13s ", corner_name );

	if( defn.Model == ModelTypeGeographic )
	{
		if (dec_flag)
		{
			printf( "(%s,", DecToDDec( x, "Long", 7 ) );
			printf( "%s)\n", DecToDDec( y, "Lat", 7 ) );
		} else {
			printf( "(%s,", GTIFDecToDMS( x, "Long", 2 ) );
			printf( "%s)\n", GTIFDecToDMS( y, "Lat", 2 ) );
		}
	} else {
		printf( "(%12.3f,%12.3f)", x, y );

		if( GTIFProj4ToLatLong( &defn, 1, &x, &y ) )
		{
			if (dec_flag)
			{
				printf( "  (%s,", DecToDDec( x, "Long", 7 ) );
				printf( "%s)", DecToDDec( y, "Lat", 7 ) );
			} else {
				printf( "  (%s,", GTIFDecToDMS( x, "Long", 2 ) );
				printf( "%s)", GTIFDecToDMS( y, "Lat", 2 ) );
			}
		}
		printf( "\n" );
	}

	if( inv_flag && GTIFPCSToImage( gtif, &x_saved, &y_saved ) ) {
		printf( "	 inverse (%11.3f,%11.3f)\n", x_saved, y_saved );
	}
	return true;
}

String GeoImage::getMouseStr(double x, double y)
{									  
	double x_saved, y_saved;
	char   str[256];
	int    n = 0;
	String retStr;

	if( !gtif )
		return false;

	x_saved = x;
	y_saved = y;								

	n += sprintf(str + n,"%s ", modelSTR.c_str() );
	n += sprintf(str + n,"(%s): ", gcsSTR.c_str() );
	if( defn.Model != ModelTypeGeographic )
	{
		n += sprintf( str + n, "[%12.7f, %12.7f]", x, y );

		if( !GTIFProj4ToLatLong( &defn, 1, &x, &y ) ) {
			return str;
		}
		n += sprintf( str + n, " || ");
	}
	n += sprintf( str + n, "[%5.7f,", x);
	n += sprintf( str + n, " %5.7f]", y);
	n += sprintf( str + n, " -- ");
	n += sprintf( str + n, "[%s,", GTIFDecToDMS( x, "Long", 2 ) );
	n += sprintf( str + n, " %s]", GTIFDecToDMS( y, "Lat", 2 ) );

	return str;
}

bool GeoImage::PrintCorners(int inv_flag, int dec_flag)
{
	printf( "Corner Coordinates:\n" );
	if( !ReportACorner( "Upper Left", 0.0, 0.0, inv_flag, dec_flag ) )
	{
		printf( " ... unable to transform points between pixel/line and PCS space\n" );
		return false;
	}

	ReportACorner( "Lower Left", 0.0, height, inv_flag, dec_flag );
	ReportACorner( "Upper Right", width, 0.0, inv_flag, dec_flag );
	ReportACorner( "Lower Right", width, height, inv_flag, dec_flag );
	ReportACorner( "Center", width/2.0, height/2.0, inv_flag, dec_flag );

	return true;

}

bool GeoImage::geoToImage(double &x, double &y)
{
	if( gtif ){ 
		// try to transform the coordinate into PCS space 
		if( !GTIFPCSToImage( gtif, &x, &y ) )
			return false;							

		if( defn.Model == ModelTypeGeographic )
			return true;
		else return GTIFProj4FromLatLong( &defn, 1, &x, &y );
	}

	x = (x - m_tieX) / m_resX + img_tieX;	
	y = (y - m_tieY) / m_resY + img_tieY;	

	return true;
}

void GeoImage::getTiePoints(double &x, double &y) 
{ 	
	x = m_tieX - img_tieX * m_resX; 
	y = m_tieY - img_tieY * m_resY; 

	if( gtif ) {
		if( defn.Model == ModelTypeGeographic )
			return;
		else 
			GTIFProj4ToLatLong( &defn, 1, &x, &y );
	}
}

bool GeoImage::imageToGeo(double &x, double &y)
{
	if( gtif ) {
		// try to transform the coordinate into PCS space 
		if( !GTIFImageToPCS( gtif, &x, &y ) )
			return false;							

		if( defn.Model == ModelTypeGeographic )
			return true;
		else return GTIFProj4ToLatLong( &defn, 1, &x, &y );
	}

	x = (x - img_tieX) * m_resX + m_tieX;
	y = (y - img_tieY) * m_resY + m_tieY;

	return true;

}

char *DecToDDec( double dfAngle, const char * pszAxis, int nPrecision )
{
	char        szFormat[30];
	static char szBuffer[50];
	const char *pszHemisphere = NULL;

	if( EQUAL(pszAxis,"Long") && dfAngle < 0.0 )
		pszHemisphere = "W";
	else if( EQUAL(pszAxis,"Long") )
		pszHemisphere = "E";
	else if( dfAngle < 0.0 )
		pszHemisphere = "S";
	else
		pszHemisphere = "N";

	sprintf( szFormat, "%%%d.%df%s", nPrecision+5, nPrecision, pszHemisphere );
	sprintf( szBuffer, szFormat, dfAngle );

	return( szBuffer );
}

bool GeoImage::readWorldFile()
{
	double rot1, rot2;
	FILE *stream;
	char *wld_filename;
	int i=0;
	char buffer[256];

	if( wf_filename == "" ) {
		wld_filename = g_strdup(filename.c_str());

		strcpy(strrchr(wld_filename, '.'), ".wld");
		stream = fopen(wld_filename, "r");
		if(!stream) {
			strcpy(strrchr(wld_filename, '.'), ".tfw");
			stream = fopen(wld_filename, "r");
			if(!stream) {
				strcpy(strrchr(wld_filename, '.'), ".jgw");
				stream = fopen(wld_filename, "r");
				if(!stream) {
					strcpy(strrchr(wld_filename, '.'), ".gfw");
					stream = fopen(wld_filename, "r");
					if(!stream) {   
						//printf("Unable to open a world file for reading\n");
						g_free(wld_filename);
						return false;
					}
				}
			}
		}

		// save world file name
		wf_filename = wld_filename;
		g_free(wld_filename);
	} else {
		stream = fopen(wf_filename.c_str(), "r");

		if(!stream) {   
			printf("Unable to open the world file for reading: %s\n", wf_filename.c_str());
			perror("fopen");
			return false;
		}
	}

	while(fgets(buffer, 255, stream)) 
	{    
		switch(i) {
			case 0:
				m_resX = atof(buffer);
				break;
			case 1:
				rot1 = atof(buffer);
				if(rot1 != 0.0)
					cout << "World file rotation 1: " << rot1 << endl;
				break;
			case 2:
				rot2 = atof(buffer);
				if(rot2 != 0.0)
					cout << "World file rotation 2: " << rot2 << endl;
				break;
			case 3:
				m_resY = atof(buffer);
				break;
			case 4:
				m_tieX = atof(buffer);
				break;
			case 5:
				m_tieY = atof(buffer);
				break;
			default:
				cout << "!!! World file has too many lines" << endl;
				fclose(stream);
				return false;
		}
		i++;
	}

	// close file and free mem
	fclose(stream);

	// set the defn.Model as geographic for point conversions
	defn.Model     = ModelTypeGeographic;

	// setup definition type for PROJ.4 conversions
	//imageProj = pj_init_plus("+proj=latlong +ellps=GRS80 +datum=GRS80");
		
	gcsSTR     = "NAD83";
	modelSTR   = "GEO-LL";

	// these images usually specific the cordinate at the center of the top-left pixel
	img_tieX   = 0.5;
	img_tieY   = 0.5;

	return TRUE;
}

void GeoImage::printImageInfo()
{
	cout << "   Image Information:" << endl;
	cout << "      " << "Model      : " << gcsSTR << endl;
	cout << "      " << "GCS        : " << modelSTR << endl;
	cout << "      " << "Dimensions : " << width << setw(10) << height << endl;
	cout << "      Tie Points --" << endl
		<< "             IMG : " << setprecision (15) << img_tieX << setw(22) << img_tieY << endl
		<< "             GEO : " << setprecision (15) << m_tieX << setw(22) << m_tieY << endl;
	cout << "      " << "Resolution --" << endl
		<< "             IMG : " << setprecision (15) << m_resX << setw(22) << m_resY << endl
		<< "            CALC : " << setprecision (15) << dx_dPp << setw(22) << dy_dPp << endl;
}
